/*!
  \file gvl_file.c
 
  \brief OGSF library - loading and manipulating volumes (lower level functions)
 
  GRASS OpenGL gsurf OGSF Library 

  (C) 1999-2008 by the GRASS Development Team
 
  This program is free software under the 
  GNU General Public License (>=v2). 
  Read the file COPYING that comes with GRASS
  for details.

  \author Tomas Paudits (February 2004)
  \author Doxygenized by Martin Landa <landa.martin gmail.com> (May 2008)
*/

#include <string.h>
#include <stdlib.h>

#include <grass/gis.h>
#include <grass/gstypes.h>
#include <grass/gsurf.h>
#include <grass/G3d.h>
#include <grass/glocale.h>

#define LUCKY                 33

#define MODE_DIRECT            0
#define MODE_SLICE             1
#define MODE_FULL              2
#define MODE_PRELOAD           3

#define MODE_DEFAULT		   0

#define STATUS_READY           0
#define STATUS_BUSY            1

/*!
  \brief structure for slice mode reading from volume file
*/
typedef struct
{
    int num, skip;
    int crnt, base;

    void *slice[MAX_VOL_SLICES];
} slice_data;

static geovol_file *Data[MAX_VOL_FILES];
static geovol_file Df[MAX_VOL_FILES];    /* trying to avoid allocation */

static int Numfiles = 0;

static int Cur_id = LUCKY;
static int Cur_max;

static int Rows, Cols, Depths;

/* local functions prototypes */
void *open_g3d_file(const char *, IFLAG *, double *, double *);
int close_g3d_file(void *);

/*!
  \brief Initialize volume files
  
  \return 1
*/
static int init_volfiles(void)
{
    int i;
    G3D_Region *w3;
    
    for (i = 0; i < MAX_VOL_FILES; i++) {
	/* avoiding dynamic allocation */
	Data[i] = &(Df[i]);
    }
    
    Cur_max = MAX_VOL_FILES;
    
    /* get window */
    w3 = GVL_get_window();
    
    /* set cols, rows, depths from window */
    Cols = w3->cols;
    Rows = w3->rows;
    Depths = w3->depths;
    
    return (1);
}

/*!
  \brief Check number of files

  \return 0
*/
static int check_num_volfiles(void)
{
    if (Numfiles < Cur_max) {
	return (0);
    }
    
    G_fatal_error (_("Maximum number of datafiles exceeded"));

    /* This return statement keeps compilers happy, it is never executed */
    return (0);
}

/*!
  \brief Get geovol_file structure for given handle

  \param id

  \return pointer to geovol_file struct
  \return NULL on failure
*/
geovol_file *gvl_file_get_volfile(int id)
{
    int i;
    
    for (i = 0; i < Numfiles; i++) {
	if (Data[i]->data_id == id) {
	    return (Data[i]);
	}
    }
    
    return (NULL);
}

/*!
  \brief Find file with name and type in geovol_file array an return handle
  
  \param name file name
  \para begin

  \param data id
  \param -1 not found
*/
int find_datah(const char *name, IFLAG type, int begin)
{
    static int i;
    int start;
    
    start = begin ? 0 : i + 1;
    
    for (i = start; i < Numfiles; i++) {
	if (!strcmp(Data[i]->file_name, name)) {
	    if (Data[i]->file_type == type) {
		return (Data[i]->data_id);
	    }
	}
    }
    
    return (-1);
}

/*!
  \brief Get file name for given handle

  \param id handle id

  \return file name
  \return NULL on failure
*/
char *gvl_file_get_name(int id)
{
    int i;
    geovol_file *fvf;
    static char retstr[GPATH_MAX];
    
    for (i = 0; i < Numfiles; i++) {
	if (Data[i]->data_id == id) {
	    fvf = Data[i];
	    strcpy(retstr, fvf->file_name);
	    
	    return (retstr);
	}
    }
    
    return (NULL);
}

/*!
  \brief Get file type for given handle

  \param vf pointer to geovol_file struct

  \return file type 
*/
int gvl_file_get_file_type(geovol_file *vf)
{
    return (vf->file_type);
}

/*!
  \brief Get data type for given handle

  \param vf pointer to geovol_file struct

  \return data type
*/
int gvl_file_get_data_type(geovol_file *vf)
{
    return (vf->data_type);
}

/*!
  \brief Get minimum and maximum value in volume file

  \param vf pointer to geovol_file struct
  \param[out] min min value
  \param[out] max max value
*/
void gvl_file_get_min_max(geovol_file *vf, double *min, double *max)
{
    *min = vf->min;
    *max = vf->max;
}

/*!
  \brief Open volume file

  \param name file name
  \param file_type file type
  \param data_type data type
  \param[out] min min value
  \param[out] max max value

  \return pointer to file
  \return NULL on failure
*/
void *open_volfile(const char *name, IFLAG file_type, IFLAG *data_type, double *min, double *max)
{
     if (file_type == VOL_FTYPE_G3D) {
        return open_g3d_file(name, data_type, min, max);
    }

    return (NULL);
}

/*!
  \brief Close volume file

  \param map volume filename
  \param type file type

  \return
  \return -1 on failure
*/
int close_volfile(void *map, IFLAG type)
{
     if (type == VOL_FTYPE_G3D) {
        return close_g3d_file(map);
    }

    return (-1);
}

/*!
  \brief Get handle for given file name and type

  \param name volume filename
  \param file_type file type

  \return data id
  \return -1 on failure
*/
int gvl_file_newh(const char *name, IFLAG file_type)
{
    geovol_file *new;
    static int first = 1;
    int i, id;
    void *m;
    IFLAG data_type;
    double min, max;

    if (first) {
	if (0 > init_volfiles()) {
	    return (-1);
	}
	
	first = 0;
    }
    
    if ( 0 <= (id = find_datah(name, file_type, 1))) {
        for (i = 0; i < Numfiles; i++) {
	    if (Data[i]->data_id == id) {
		new = Data[i];
		new->count++;
		
		return (id);
	    }
        }
    }
    
    if (0 > check_num_volfiles()) {
	return (-1);
    }

    if (!name) {
	return (-1);
    }
    
    if ((m = open_volfile(name, file_type, &data_type, &min, &max)) == NULL) {
	return (-1);
    }
    
    new = Data[Numfiles];
    
    if (new) {
	Numfiles++;
	new->data_id = Cur_id++;
	
	new->file_name = G_store(name);
	new->file_type = file_type;
	new->count = 1;
	new->map = m;
	new->min = min;
	new->max = max;
	new->data_type = data_type;
	
	new->status = STATUS_READY;
	new->buff = NULL;
	
	new->mode = 255;
	gvl_file_set_mode(new, MODE_DEFAULT);
	
	return (new->data_id);
    }
    
    return (-1);
}

/*!
  \brief Free allocated buffers

  \param vf pointer to geovol_file struct

  \return 1
*/
int free_volfile_buffs(geovol_file * vf)
{
    if (vf->mode == MODE_SLICE) {
        G_free(vf->buff);
        vf->buff = NULL;
    }
    
    if (vf->mode == MODE_PRELOAD) {
	G_free(vf->buff);
	vf->buff = NULL;
    }
    
    return (1);
}

/*!
  \brief Free geovol_file structure for given handle

  \param id 

  \return
*/
int gvl_file_free_datah(int id)
{
    int i, j, found = -1;
    geovol_file *fvf;

    G_debug (3, "gvl_file_free_datah");

    for (i = 0; i < Numfiles; i++) {
	if (Data[i]->data_id == id) {
	    found = 1;
	    fvf = Data[i];
	    
	    if (fvf->count > 1) {
		fvf->count--;
	    } else {
		close_volfile(fvf->map, fvf->file_type);
		free_volfile_buffs(fvf);

		G_free(fvf->file_name);
		fvf->file_name = NULL;

		fvf->data_id = 0;
		
		for (j = i; j < (Numfiles - 1); j++) {
		    Data[j] = Data[j + 1];
		}
		
		Data[j] = fvf;
		
		--Numfiles;
	    }
	}
    }
    
    return (found);
}

/******************************************************************/
/* reading from G3D raster volume files */
/******************************************************************/

/*!
  \brief Open g3d file

  \param filename file name
  \param type data type
  \param[out] min min value
  \param[out] max max value

  \pointer to data
*/
void *open_g3d_file(const char *filename, IFLAG *type, double *min, double *max)
{
    const char *mapset;
    int itype;
    void *map;

    /* search for g3d file a return his mapset */
    mapset = G_find_grid3(filename, "");
    if (!mapset) {
	G_warning (_("3D raster map <%s> not found"),
		   filename);
        return (NULL);
    }

    /* open g3d file */
    map = G3d_openCellOld(filename, mapset, G3D_DEFAULT_WINDOW, G3D_TILE_SAME_AS_FILE, G3D_USE_CACHE_DEFAULT);
    if (!map) {
        G_warning (_("Unable to open 3D raster map <%s>"), filename);
        return (NULL);
    }

    /* load range into range structure of map */
    if (!G3d_range_load(map)) {
        G_warning ("Unable to read range of 3D raster map <%s>", filename);
        return (NULL);
    }

    G3d_range_min_max(map, min, max);
    
    /* get file data type */
    itype = G3d_fileTypeMap(map);
    if (itype == FCELL_TYPE)
        *type = VOL_DTYPE_FLOAT;
    if (itype == DCELL_TYPE)
        *type = VOL_DTYPE_DOUBLE;
    
    return (map);
}

/*!
  \brief Close g3d file

  \param map 3d raster map

  \return -1 on failure
  \return 1 on success
*/
int close_g3d_file(void *map)
{
    /* close opened g3d file */
    if (G3d_closeCell((G3D_Map *)map) != 1) {
	G_warning (_("Unable to close 3D raster map <%s>"), ((G3D_Map *)map)->fileName);
	return (-1);
    }
    
    return (1);
}

/*!
  \brief Eead value from g3d file

  \param type data type
  \param map 3D raster map
  \param x,y,z real coordinates
  \param[out] value data value

  \return -1 on failure
  \return 1 on success
*/
int read_g3d_value(IFLAG type, void *map, int x, int y, int z, void *value)
{
    switch (type) {
        /* float data type */
    case (VOL_DTYPE_FLOAT) :
	*((float *)value) = G3d_getFloat(map, x, y, z);
	break;
	
        /* double data type */
    case (VOL_DTYPE_DOUBLE) :
	*((double *)value) = G3d_getDouble(map, x, y, z);
	break;
	
        /* unsupported data type */
    default :
	return (-1);
    }
    
    return (1);
}

/*!
  \brief Read slice of values at level from g3d file

  \param type data type
  \param map 3D raster map
  \param level
  \param[out] data

  \return -1 on failure
  \return 0 on success
*/
int read_g3d_slice(IFLAG type, void *map, int level, void *data)
{
    int x, y;
    
    switch (type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            for (x = 0; x < Cols; x++) {
                for (y = 0; y < Rows; y++) {       
                    ((float *) data)[x + y * Cols] =
                        G3d_getFloat(map, x, y, level);
                }
            }       
        
            break;
        
        /* double data type */
        case (VOL_DTYPE_DOUBLE) :                        
            for (x = 0; x < Cols; x++) {
                for (y = 0; y < Rows; y++) {       
                    ((double *) data)[x + y * Cols] =
                        G3d_getDouble(map, x, y, level);
                }
            }       
            
            break; 
        
        /* unsupported data type */
        default :
            return (-1);
    }        
        
    return (1);
}

/*!
  \brief Read all values from g3d file

  \param type data type
  \param map 3D raster map
  \param[out] data data buffer

  \return -1 on failure
  \return 1 on success
*/
int read_g3d_vol(IFLAG type, void *map, void *data)
{    
    int x, y, z;
    
    switch (type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            for (x = 0; x < Cols; x++) {
                for (y = 0; y < Rows; y++) {
                    for (z = 0; z < Depths; z++) {
                        ((float *) data)[x + y * Cols + z * Rows * Cols] =
                            G3d_getFloat(map, x, y, z);
                    }
                }
            }       
        
            break;
        
        /* double data type */
        case (VOL_DTYPE_DOUBLE) :                        
            for (x = 0; x < Cols; x++) {
                for (y = 0; y < Rows; y++) {
                    for (z = 0; z < Depths; z++) {
                        ((double *) data)[x + y * Cols + z * Rows * Cols] =
                            G3d_getDouble(map, x, y, z);
                    }
                }
            }       
            
            break; 
        
        /* unsupported data type */
        default :
            return (-1);
    }        
        
    return (1);
}

/*!
  \brief Check for null value

  \param type data type
  \param value

  \return 1 if value is null
  \return 0 if value is not null
  \return -1 on failure (unsupported data type
*/
int is_null_g3d_value(IFLAG type, void *value)
{
    switch (type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            return G3d_isNullValueNum(value, FCELL_TYPE);
            break;

        /* double data type */
        case (VOL_DTYPE_DOUBLE) :
            return G3d_isNullValueNum(value, DCELL_TYPE);
            break;

        /* unsupported data type */
        default :
            return (-1);
    }

    return (-1);
}

/******************************************************************/
/* reading from buffer */
/******************************************************************/

/*!
  \brief Get value from buffer

  \param type data type
  \param data data buffer
  \param offset
  \param value

  \return -1 on failure (unsupported data type)
  \return 1 on success
*/
int get_buff_value(IFLAG type, void *data, int offset, void *value)
{
    switch (type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            *((float*) value) = ((float *) data)[offset];
            break;

        /* double data type */
        case (VOL_DTYPE_DOUBLE) :
            *((double*) value) = ((double *) data)[offset];
            break;

        /* unsupported data type */
        default :
            return (-1);
    }

    return (1);
}

/******************************************************************/
/* direct mode reading from volume file */
/******************************************************************/

/*!
  \brief Read value direct from volume file

  \param vf pointer to geovol_file struct
  \param x,y,z real point
  \oaram[out] value data value

  \return -1 on failure
  \return 1 on success
*/
int get_direct_value(geovol_file *vf, int x, int y, int z, void *value)
{
    switch (vf->file_type) {
        /* G3D file type */
        case (VOL_FTYPE_G3D) :
            if (0 > read_g3d_value(vf->data_type, vf->map, x, y, z, value))
                return (-1);
            break;

        default :
            return (-1);
    }

    return (1);
}

/******************************************************************/
/* full mode reading from volume file */
/******************************************************************/

/*!
  \brief Allocate buffer memory for full mode reading

  \param vf pointer to geovol_file

  \return -1 on failure
  \return 1 on success
*/
int alloc_vol_buff(geovol_file *vf)
{
    switch (vf->data_type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            if ((vf->buff = (float*) G_malloc(sizeof(float) * Cols * Rows * Depths)) == NULL)
                return (-1);
            break;

        /* double data type */
        case (VOL_DTYPE_DOUBLE) :
            if ((vf->buff = (double*) G_malloc(sizeof(double) * Cols * Rows * Depths)) == NULL)
                return (-1);
            break;

        /* unsupported data type */
        default :
            return (-1);
    }

    return (1);
}

/*!
  \brief Free memory buffer memory

  \param vf pointer to geovol_file struct

  \return 1
*/
int free_vol_buff(geovol_file *vf)
{
    G_free(vf->buff);

    return (1);
}

/*!
  \brief Read all values from volume file

  \param vf pointer to geovol_file struct

  \return -1 on failure
  \return 1 on success
*/
int read_vol(geovol_file *vf)
{
    switch (vf->file_type) {
        /* G3D file format */
        case (VOL_FTYPE_G3D) :
            if (0 > read_g3d_vol(vf->data_type, vf->map, vf->buff))
                    return (-1);
            break;
        /* unsupported file format */
        default :
            return (-1);
    }

    return (1);
}

/*!
  \brief Get value from volume buffer

  \param vf pointer to geovol_file struct
  \param x,y,z real point
  \param[out] value data value

  \return 1
*/
int get_vol_value(geovol_file *vf, int x, int y, int z, void *value)
{
    get_buff_value(vf->data_type, vf->buff, z * Rows * Cols + y * Cols + x, value);

    return (1);
}

/******************************************************************/
/* slice mode reading from volume file */
/******************************************************************/

/*!
  \brief Allocate buffer for slice mode reading
  
  \param vf pointer to geovol_file struct

  \return -1 on failure
  \return 1 on success
*/
int alloc_slice_buff(geovol_file *vf)
{
    int i;
    slice_data *sd = (slice_data *) vf->buff;

    switch (vf->data_type) {
        /* float data type */
        case (VOL_DTYPE_FLOAT) :
            for (i = 0; i < sd->num; i++) {
                if ((sd->slice[i] = (float*) G_malloc(sizeof(float) * Cols * Rows)) == NULL)
                    return (-1);
            }
            break;

        /* double data type */
        case (VOL_DTYPE_DOUBLE) :
            for (i = 0; i < sd->num; i++) {
                if ((sd->slice[i] = (double*) G_malloc(sizeof(double) * Cols * Rows)) == NULL)
                    return (-1);
            }
            break;

        /* unsupported data type */
        default :
            return (-1);
    }

    return (1);
}

/*!
  \brief Free buffer for slice mode reading

  \param vf pointer to geovol_file struct

  \return 1
*/
int free_slice_buff(geovol_file *vf)
{
    int i;
    slice_data *sd = (slice_data *) vf->buff;

    for (i = 0; i < sd->num; i++) {
        G_free(sd->slice[i]);
    }

    return (1);
}

/*!
  \brief Read slice of values at level from volume file

  \param vf pointer to geovol_file struct
  \param s
  \param l

  \return -1 on failure
  \return 1 on success
*/
int read_slice(geovol_file *vf, int s, int l)
{
    /* get slice structure */
    slice_data *sd = (slice_data *) vf->buff;

    switch (vf->file_type) {
    /* G3D file format */
        case (VOL_FTYPE_G3D) :
            if (0 > read_g3d_slice(vf->data_type, vf->map, l, sd->slice[s]))
		return (-1);
            break;
        /* unsupported file format */
        default :
            return (-1);
    }

    return (1);
}

/*!
  \brief Read new slice into buffer

  \param vf pointer to geovol_file struct
*/
void shift_slices(geovol_file *vf)
{
    void *tmp;
    int i;

    slice_data *sd = (slice_data *) vf->buff;

    /* change pointers to slices */
    tmp = sd->slice[0];
    for (i = 0; i < sd->num - 1; i++) {
        sd->slice[i] = sd->slice[i+1];
    }
    sd->slice[sd->num - 1] = tmp;

    /* read new slice data */
    read_slice(vf, sd->num, sd->crnt + 1 + (sd->num - sd->base));

    /* increase current slice value */
    sd->crnt++;
}

/*!
  \brief Get value from slice buffer

  \param vf pointer to geovol_file struct
  \param x,y,z real point
  \param[out] value data value

  \return -1 on failure
  \return 1 on success
*/
int get_slice_value(geovol_file *vf, int x, int y, int z, void *value)
{
    slice_data *sd = (slice_data *) vf->buff;

    /* value is in loaded slices */
    if ((z >= sd->crnt - (sd->base - 1)) &&
        (z <= sd->crnt + sd->num - sd->base)) {

        get_buff_value(vf->data_type, sd->slice[(z - sd->crnt)], y * Cols + x, value);
    }

    /* if value isn't in loaded slices, need read new data slice */
    else if (z == sd->crnt - (sd->base - 1) + 1) {
        shift_slices(vf);
        get_buff_value(vf->data_type, sd->slice[(z - sd->crnt)], y * Cols + x, value);
    }

    /* value out of range */
    else {
        return (-1);
    }

    return (1);
}

/******************************************************************/
/* reading from volume file */
/******************************************************************/

/*!
  \brief Start read - allocate memory buffer a read first data into buffer

  \param vf pointer to geovol_file struct

  \return -1 on failure
  \return 1 on success
*/
int gvl_file_start_read(geovol_file *vf)
{
    int i;
    slice_data *sd;

    /* check status */
    if (vf->status == STATUS_BUSY)
         return (-1);
         
    switch (vf->mode) {
        /* read whole volume into memory */
        case (MODE_FULL) :
             /* allocate memory */
            if (0 > alloc_vol_buff(vf))
                return (-1);

            /* read volume */
            read_vol(vf);
            break;

        /* read first data slices into memory */
        case (MODE_SLICE) :
            /* allocate slices buffer memory */
            if (0 > alloc_slice_buff(vf))
                return (-1);

            /* read volume */
            sd = (slice_data *) vf->buff;
            /* set current slice to 0 */
            sd->crnt = 0;

            /* read first slices into buffer */
            for (i = 0; i < (sd->num - sd->base + 1); i++)
                 read_slice(vf, (sd->base - 1) + i, i);
            break;
     }

     /* set status */
     vf->status = STATUS_BUSY;

    return (1);
}

/*!
  \brief End read - free buffer memory

  \param vf pointer to geovol_file struct

  \return -1 on failure
  \return 1 on success
*/
int gvl_file_end_read(geovol_file *vf)
{
     /* check status */
     if (vf->status == STATUS_READY)
         return (-1);

    switch (vf->mode) {
        case (MODE_FULL) :
             if (0 > free_vol_buff(vf))
                return (-1);
            break;

        case (MODE_SLICE) :
            /* allocate slices buffer memory */
            if (0 > free_slice_buff(vf))
                return (-1);
     }

     /* set status */
     vf->status = STATUS_READY;

    return (1);
}

/*!
  \brief Get value for volume file at x, y, z

  \param vf pointer to geovol_file struct

  \return -1 on failure
  \return 1 on success
*/
int gvl_file_get_value(geovol_file *vf, int x, int y, int z, void *value)
{
    /*  check status */
    if (vf->status != STATUS_BUSY) {
        return (-1);
    }

    switch (vf->mode) {
        /* read value direct from g3d file */
        case (MODE_DIRECT) :
            if (0 > get_direct_value(vf, x, y, z, value))
                return (-1);
            break;

        case (MODE_SLICE) :
            if (0 > get_slice_value(vf, x, y, z, value))
                return (-1);
            break;

        case (MODE_FULL) :
		case (MODE_PRELOAD) :
            if (0 > get_vol_value(vf, x, y, z, value))
                return (-1);
            break;
    }
    return (1);
}

/*!
  \brief Check for null value

  \param vf pointer to geovol_file struct
  \param value data value

  \return -1 on failure
  \return 1 on success
*/  
int gvl_file_is_null_value(geovol_file *vf, void *value)
{
    switch (vf->file_type) {
        /* G3D file format */
        case (VOL_FTYPE_G3D) :
            return is_null_g3d_value(vf->file_type, value);
            break;
        /* unsupported file format */
        default :
            return (-1);
    }

    return (-1);
}

/******************************************************************/
/* set parameters for reading volumes */
/******************************************************************/

/*!
  \brief Set read mode

  \param vf pointer to geovol_file struct
  \param mode

  \return -1 on failure
  \return 1 on success
*/
int gvl_file_set_mode(geovol_file *vf, IFLAG mode)
{
    slice_data *sd;
    
    if (vf->status == STATUS_BUSY)
        return (-1);
    
    if (vf->mode == mode)
        return (1);
    
    if (vf->mode == MODE_SLICE)
        G_free(vf->buff);
    
    if (vf->mode == MODE_PRELOAD)
	G_free(vf->buff);
    
    if (mode == MODE_SLICE) {
	if ((vf->buff = (slice_data *) G_malloc(sizeof(slice_data))) == NULL)
		return (-1);
	
	sd = (slice_data *) vf->buff;
	sd->num = 1;
	sd->crnt = 0;
	sd->base = 1;
    }
    
    if (mode == MODE_PRELOAD) {
	/* allocate memory */
	if (0 > alloc_vol_buff(vf))
	    return (-1);
	
	/* read volume */
	read_vol(vf);
    }
    
    vf->mode = mode;

    return (1);
}

/*!
  \brief Set parameters for slice reading

  \param vf pointer to geovol_file struct
  \param n
  \param b

  \return -1 on failure
  \return 1 on success
*/
int gvl_file_set_slices_param(geovol_file *vf, int n, int b)
{
    slice_data *sd;

    if (vf->status == STATUS_BUSY)
        return (-1);

    if (!(vf->mode == MODE_SLICE))
        return (-1);

    sd = (slice_data *) vf->buff;
    sd->num = n;
    sd->base = b;

    return (1);
}
